!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine write_vtk_rectilinear (datavtk_name, var_name, var, x, y, z, nx, ny, nz, n_var, var_dim, format, header)
!==================================================================================================================================|

!----------------------------------------------------------------------------------------------------------------------------------|
!  Example :
!  datavtk_name = 'output_filename.vtk' (character, scalar, the output filename)
!  var_name = (/ 'density', 'velocity', 'pressure' /) (character, array, the name of variables)
!  var(1, :, :, :) = ro (real*4, array, the density data)
!  var(2, :, :, :) = vx (real*4, array, the x-component of velocity data)
!  var(3, :, :, :) = vy (real*4, array, the y-component of velocity data)
!  var(4, :, :, :) = vz (real*4, array, the z-component of velocity data)
!  var(5, :, :, :) = pr (real*4, array, the pressure data)
!  x (real*8, array, the x-coordinate)
!  y (real*8, array, the y-coordinate)
!  z (real*8, array, the z-coordinate)
!  nx (integer*4, scalar, the number of x-coordinate points)
!  ny (integer*4, scalar, the number of y-coordinate points)
!  nz (integer*4, scalar, the number of z-coordinate points)
!  n_var = 3 (integer*4, scalar, the number of variables)
!  var_dim = (/1, 3, 1/) (integer*4, array, the number of components of variables)
!  formate = 0 (integer*4, scalar, the output format, 0 is unformatted, 1 is formatted)
!  header = 'Really cool data' (character, scalar, the information of your data, up to 256 characters) 
!----------------------------------------------------------------------------------------------------------------------------------|

   implicit none

   integer(4), parameter :: n_dim = 3
   integer(4), parameter :: n_point_cell = 8
   integer(4), parameter :: file_unit = 20

   integer(4), intent(in) :: n_var, format
   integer(4), intent(in) :: nx, ny, nz

   character(*), intent(in) :: var_name(n_var), datavtk_name
   character(256), intent(in) :: header

   integer(4), intent(in) :: var_dim(n_var)
   real(8), intent(in) :: x(nx), y(ny), z(nz)
   real(4), intent(in) :: var(sum(var_dim), nx, ny, nz)

   integer(4) :: i_var, index_var

!----------------------------------------------------------------------------------------------------------------------------------|
!  write the head data
!----------------------------------------------------------------------------------------------------------------------------------|
   open (file_unit, file = trim (datavtk_name), form = 'formatted')
   write (file_unit, '(a)') '# vtk DataFile Version 3.0'
   write (file_unit, '(a)') trim(header)

   if (format .eq. 0) then
      write (file_unit, '(a)') 'BINARY'
   else
      write (file_unit, '(a)') 'ASCII'
   endif
   close (file_unit)

!----------------------------------------------------------------------------------------------------------------------------------|
!  write coordinates
!----------------------------------------------------------------------------------------------------------------------------------|
   open (file_unit, file = trim (datavtk_name), form = 'formatted', position = 'append')
   write (file_unit, '(a)') 'DATASET RECTILINEAR_GRID'
   write (file_unit, '(a, 1x, i, 1x, i, 1x, i)') 'DIMENSIONS', nx, ny, nz
   write (file_unit, '(a, 1x, i, 1x, a)') 'X_COORDINATES', nx, 'double'
   close (file_unit)

   if (format .eq. 0) then
      open (file_unit, file = trim (datavtk_name), form = 'unformatted', position = 'append',                                      &
         convert = 'big_endian', access = 'stream')
      write (file_unit) x
      close (file_unit)
   else
      open (file_unit, file = trim (datavtk_name), form = 'formatted', position = 'append')
      write (file_unit, '(6es)') x
      close (file_unit)
   endif

   open (file_unit, file = trim (datavtk_name), form = 'formatted', position = 'append')
   write (file_unit, '(a, 1x, i, 1x, a)') 'Y_COORDINATES', ny, 'double'
   close (file_unit)

   if (format .eq. 0) then
      open (file_unit, file = trim (datavtk_name), form = 'unformatted', position = 'append',                                      &
         convert = 'big_endian', access = 'stream')
      write (file_unit) y
      close (file_unit)
   else
      open (file_unit, file = trim (datavtk_name), form = 'formatted', position = 'append')
      write (file_unit, '(6es)') y
      close (file_unit)
   endif

   open (file_unit, file = trim (datavtk_name), form = 'formatted', position = 'append')
   write (file_unit, '(a, 1x, i, 1x, a)') 'Z_COORDINATES', ny, 'double'
   close (file_unit)

   if (format .eq. 0) then
      open (file_unit, file = trim (datavtk_name), form = 'unformatted', position = 'append',                                      &
         convert = 'big_endian', access = 'stream')
      write (file_unit) z
      close (file_unit)
   else
      open (file_unit, file = trim (datavtk_name), form = 'formatted', position = 'append')
      write (file_unit, '(6es)') z
      close (file_unit)
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  write the point data
!----------------------------------------------------------------------------------------------------------------------------------|
   open (file_unit, file = trim (datavtk_name), form = 'formatted', position = 'append')
   write (file_unit, '(a, i)') 'POINT_DATA', nx * ny * nz
   close (file_unit)

   index_var = 1
   do i_var = 1, n_var
!----------------------------------------------------------------------------------------------------------------------------------|
!     write the scalar point data
!----------------------------------------------------------------------------------------------------------------------------------|
      if (var_dim(i_var) .eq. 1) then
         open (file_unit, file = trim (datavtk_name), form = 'formatted', position = 'append')
         write (file_unit, '(a, 1x, a, 1x, a)') 'SCALARS', trim (var_name(i_var)), 'float'
         write (file_unit, '(a)') 'LOOKUP_TABLE default'
         close (file_unit)
         if (format .eq. 0) then
            open (file_unit, file = trim (datavtk_name), form = 'unformatted', position = 'append',                                &
               convert = 'big_endian', access = 'stream')
            write (file_unit) var(index_var, :, :, :)
            close (file_unit)
         else
            open (file_unit, file = trim (datavtk_name), form = 'formatted', position = 'append')
            write (file_unit, '(6es)') var(index_var, :, :, :)
            close (file_unit)
         endif
         index_var = index_var + var_dim(i_var)
      elseif (var_dim(i_var) .eq. 3) then
!----------------------------------------------------------------------------------------------------------------------------------|
!     write the vector cell data
!----------------------------------------------------------------------------------------------------------------------------------|
         open (file_unit, file = trim (datavtk_name), form = 'formatted', position = 'append')
         write (file_unit, '(a, 1x, a, 1x, a)') 'VECTORS', trim (var_name(i_var)), 'float'
         close (file_unit)
         if (format .eq. 0) then
            open (file_unit, file = trim (datavtk_name), form = 'unformatted', position = 'append',                                &
               convert='big_endian', access = 'stream')
            write (file_unit) var(index_var : index_var + 2, :, :, :)
            close (file_unit)
         else
            open (file_unit, file = trim (datavtk_name), form = 'formatted', position = 'append')
            write (file_unit, '(6es)') var(index_var : index_var + 2, :, :, :)
            close (file_unit)
         endif
         index_var = index_var + var_dim(i_var)
      endif
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine write_vtk_rectilinear
